library(MODIStsp)
library(sf)
library(tidyverse)
library(tmap)
library(tidycensus)
library(raster)
library(viridis)
In order to use the MODIStsp package, you must register on the EARTHDATA page, at this link.
We will extract the MODIS land cover raster (MCD12C1) and make a simple map of land cover in the state of Colorado. A description of the data can be found here.
First, let’s extract a polygon of Colorado’s state borders using tidycensus.
co_polygon<-get_decennial(geography = "state",
variables = "P001001",
year = 2010,
geometry=TRUE) %>%
filter(NAME=="Colorado") %>%
st_transform(4326)
In case we want to overlay county-level boundaries against the land cover raster, we can also extract a spatial dataset of Colorado counties:
co_counties_polygon<-get_decennial(geography = "county",
state="CO",
variables = "P001001",
year = 2010,
geometry=TRUE) %>%
st_transform(4326)
After extracting this data from tidycensus, we’ll need to extract its bounding box, or alternatively, download the data to our local machine; this will allow us to specify the desired extent of the raster when we extract it using the MODIStsp function. In this case, we’ll go with the first option.
First, we’ll define a path and name for the spatial dataset that is to be downloaded:
filepath<-"land_cover/colorado.shp"
Then, we’ll download the co_polygon spatial data object as a shapefile using the st_write function:
st_write(co_polygon, paste0(filepath))
The Colorado state polygon is stored in a sub-directory (named “land_cover”) within our working directory, and is named “colorado.shp”.
Now, we’re ready to
MODIStsp_get_prodlayers("MCD12Q1")
## $prodname
## [1] "LandCover_Type_Yearly_500m (MCD12Q1)"
##
## $bandnames
## [1] "LC1" "LC2" "LC3"
## [4] "LC4" "LC5" "LC_Prop1"
## [7] "LC_Prop2" "LC_Prop3" "LC_Prop1_Assessment"
## [10] "LC_Prop2_Assessment" "LC_Prop3_Assessment" "LC_QC"
## [13] "LC_LW"
##
## $bandfullnames
## [1] "Land Cover Type 1 (IGBP)*"
## [2] "Land Cover Type 2 (UMD)*"
## [3] "Land Cover Type 3 (LAI/fPAR)*"
## [4] "Land Cover Type 4 (NPP/BGC)*"
## [5] "Land Cover Type 5: Annual Plant Functional Types classification"
## [6] "FAO-Land Cover Classification System 1 (LCCS1) land cover layer"
## [7] "FAO-LCCS2 land use layer"
## [8] "FAO-LCCS3 surface hydrology layer"
## [9] "LCCS1 land cover layer confidence"
## [10] "LCCS2 land use layer confidence"
## [11] "LCCS3 surface hydrology layer confidence"
## [12] "Land Cover QC"
## [13] "Binary land/water mask derived from MOD44W"
##
## $quality_bandnames
## NULL
##
## $quality_fullnames
## NULL
##
## $indexes_bandnames
## NULL
##
## $indexes_fullnames
## NULL
MODIStsp(gui = FALSE,
out_folder = "land_cover",
out_folder_mod = "land_cover",
selprod = "LandCover_Type_Yearly_500m (MCD12Q1)",
bandsel = "LC1",
user = "YOUR_USERNAME" ,
password = "YOUR_PASSWORD",
start_date = "2019.01.01",
end_date = "2019.12.31",
output_proj = 4326,
verbose = FALSE,
spatmeth = "file",
spafile = filepath,
out_format = "GTiff")
CO_raster<-raster("land_cover/colorado/LandCover_Type_Yearly_500m_v6/LC1/MCD12Q1_LC1_2019_001.tif")
raster_label<-c("Evergreen needleleaf forests",
"Evergreen broadleaf forests",
"Deciduous needleleaf forests",
"Deciduous broadleaf forests",
"Mixed forests",
"Closed shrublands",
"Open shrublands",
"Woody savannas",
"Savannas",
"Grasslands",
"Permanent wetlands",
"Croplands",
"Urban and built-up lands",
"Cropland/natural vegetation mosaics",
"Snow and ice",
"Barren",
"Water bodies")
CO_landcover_tmap<-tm_shape(CO_raster)+
tm_raster(breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18),
labels=raster_label,
palette="Set1",
title="Colorado Land Cover, 2019")+
tm_legend(outside=T)
tmap_mode("plot")
## tmap mode set to plotting
CO_landcover_tmap
## Some legend labels were too wide. These labels have been resized to 0.58. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
tmap_mode("view")
## tmap mode set to interactive viewing
CO_landcover_tmap
tmap_mode("plot")
## tmap mode set to plotting
CO_landcover_tmap+
tm_shape(co_counties_polygon)+
tm_polygons(alpha=0)
## Some legend labels were too wide. These labels have been resized to 0.58. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
tmap_mode("view")
## tmap mode set to interactive viewing
CO_landcover_tmap+
tm_shape(co_counties_polygon)+
tm_polygons(alpha=0)
CO_raster_df <- as.data.frame(CO_raster, xy = TRUE) %>%
mutate(MCD12Q1_LC1_2019_001 = as.factor(MCD12Q1_LC1_2019_001))
levels(CO_raster_df$MCD12Q1_LC1_2019_001)
## [1] "1" "4" "5" "7" "8" "9" "10" "11" "12" "13" "14" "15" "16" "17"
CO_raster_df$MCD12Q1_LC1_2019_001<-factor(CO_raster_df$MCD12Q1_LC1_2019_001, levels=c("1", "2",
"3", "4", "5",
"6", "7", "8",
"9", "10", "11",
"12", "13", "14",
"15", "16", "17"))
levels(CO_raster_df$MCD12Q1_LC1_2019_001) <- raster_label
levels(CO_raster_df$MCD12Q1_LC1_2019_001)
## [1] "Evergreen needleleaf forests" "Evergreen broadleaf forests"
## [3] "Deciduous needleleaf forests" "Deciduous broadleaf forests"
## [5] "Mixed forests" "Closed shrublands"
## [7] "Open shrublands" "Woody savannas"
## [9] "Savannas" "Grasslands"
## [11] "Permanent wetlands" "Croplands"
## [13] "Urban and built-up lands" "Cropland/natural vegetation mosaics"
## [15] "Snow and ice" "Barren"
## [17] "Water bodies"
ggplot() +
geom_raster(data = CO_raster_df,
aes(x = x, y = y, fill = MCD12Q1_LC1_2019_001)) +
geom_sf(data = co_polygon, inherit.aes = FALSE, fill = NA) +
scale_fill_viridis(name = "Land Cover Type", discrete = TRUE) +
labs(title = "Land Cover classification in Colorado",
subtitle = "01-01-2019 - 31-12-2019",
x = "Longitude",
y = "Latitude") +
theme_minimal()
ggplot() +
geom_raster(data = CO_raster_df,
aes(x = x, y = y, fill = MCD12Q1_LC1_2019_001)) +
geom_sf(data=co_counties_polygon, inherit.aes=F, fill=NA)+
scale_fill_viridis_d(name = "Land Cover Type", option="C", direction=-1 ) +
labs(title = "Land Cover classification in Colorado",
subtitle = "01-01-2019 - 31-12-2019",
x = "Longitude",
y = "Latitude") +
theme_minimal()
CO_landcover_ggplot<-
ggplot() +
geom_raster(data = CO_raster_df,
aes(x = x, y = y, fill = MCD12Q1_LC1_2019_001)) +
geom_sf(data=co_counties_polygon, inherit.aes=F, fill=NA)+
labs(title = "Land Cover classification in Colorado",
subtitle = "01-01-2019 - 31-12-2019",
x = "Longitude",
y = "Latitude") +
theme_minimal()
colors<-c("lightpink1", "lightgrey", "magenta", "lightcyan",
"yellow", "palegreen", "seashell", "plum",
"skyblue", "salmon2", "khaki", "tomato",
"gray0", "springgreen", "purple", "orange",
"mediumblue", "slategray")
CO_landcover_ggplot +
scale_fill_manual(values=colors, drop=F)
MODIStsp_get_prodlayers("M*D13A2")
## $prodname
## [1] "Vegetation_Indexes_16Days_1Km (M*D13A2)"
##
## $bandnames
## [1] "NDVI" "EVI" "VI_QA" "b1_Red" "b2_NIR" "b3_Blue"
## [7] "b7_SWIR" "View_Zen" "Sun_Zen" "Rel_Az" "DOY" "Rely"
##
## $bandfullnames
## [1] "16 day NDVI average" "16 day EVI average"
## [3] "VI quality indicators" "Surface Reflectance Band 1"
## [5] "Surface Reflectance Band 2" "Surface Reflectance Band 3"
## [7] "Surface Reflectance Band 7" "View zenith angle of VI pixel"
## [9] "Sun zenith angle of VI pixel" "Relative azimuth angle of VI pixel"
## [11] "Day of year of VI pixel" "Quality reliability of VI pixel"
##
## $quality_bandnames
## [1] "QA_qual" "QA_usef" "QA_aer" "QA_adj_cld" "QA_BRDF"
## [6] "QA_mix_cld" "QA_land_wat" "QA_snow_ice" "QA_shd"
##
## $quality_fullnames
## [1] "VI Quality"
## [2] "VI usefulness"
## [3] "Aerosol quantity"
## [4] "Adjacent cloud detected"
## [5] "Atmosphere BRDF correction performed"
## [6] "Mixed Clouds"
## [7] "Land/Water Flag"
## [8] "Possible snow/ice"
## [9] "Possible shadow"
##
## $indexes_bandnames
## [1] "SR" "NDFI" "NDII7" "SAVI"
##
## $indexes_fullnames
## [1] "Simple Ratio (NIR/Red)"
## [2] "Flood Index (Red-SWIR2)/(Red+SWIR2)"
## [3] "NDII7 (NIR-SWIR2)/(NIR+SWIR2)"
## [4] "SAVI (NIR-Red)/(NIR+Red+0.5)*(1+0.5)"
MODIStsp(
gui = FALSE,
out_folder = "VegetationData",
out_folder_mod = "VegetationData",
selprod = "Vegetation_Indexes_16Days_1Km (M*D13A2)",
bandsel = "NDVI",
user = "YOUR_USERNAME",
password = "YOUR_PASSWORD",
start_date = "2020.07.03",
end_date = "2020.07.03",
output_proj = 4326,
verbose = FALSE,
spatmeth = "file",
spafile = "VegetationData/colorado.shp",
out_format = "GTiff"
)
CO_ndvi<-raster("VegetationData/colorado/VI_16Days_1Km_v6/NDVI/MYD13A2_NDVI_2020_185.tif")
CO_ndvi
## class : RasterLayer
## dimensions : 269, 471, 126699 (nrow, ncol, ncell)
## resolution : 0.01490176, 0.01491085 (x, y)
## extent : -109.0603, -102.0415, 36.99243, 41.00344 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : /Users/adra7980/Documents/git_repositories/ARSC5040_GIS/general_resources/VegetationData/colorado/VI_16Days_1Km_v6/NDVI/MYD13A2_NDVI_2020_185.tif
## names : MYD13A2_NDVI_2020_185
## values : -32768, 32767 (min, max)
See here for scale factor.
CO_ndvi_adjusted<-CO_ndvi*0.0001
tmap_mode("plot")
## tmap mode set to plotting
CO_ndvi_map<-tm_shape(CO_ndvi_adjusted)+
tm_raster(palette = "YlGn", style="cont", breaks=c(-0.2, 0,0.2,0.4,0.6,0.8, 1))+
tm_legend(outside=T)
CO_ndvi_map
## Warning: Breaks contains positive and negative values. Better is to use
## diverging scale instead, or set auto.palette.mapping to FALSE.
tmap_mode("view")
## tmap mode set to interactive viewing
CO_ndvi_map
## Warning: Breaks contains positive and negative values. Better is to use
## diverging scale instead, or set auto.palette.mapping to FALSE.
NDVI_co_df <- as.data.frame(CO_ndvi_adjusted, xy = TRUE, na.rm = TRUE)
# Visualising using ggplot2
ggplot() +
geom_raster(
data = NDVI_co_df,
aes(x = x, y = y, fill = MYD13A2_NDVI_2020_185)
) +
geom_sf(data = co_counties_polygon, inherit.aes = FALSE, fill = NA) +
scale_fill_viridis(name = "NDVI") +
labs(
title = "NDVI (Normalized Difference Vegetation Index) in Colorado",
subtitle = "03-07-2020",
x = "Longitude",
y = "Latitude"
) +
theme_minimal()